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Abstract 



We perform extensive simulations of the sandpile model on a Sierpinski gasket. Critical 
exponents for waves and avalanches are determined. We extend the existing theory of waves 
to the present case. This leads to an exact value for the exponent t w which describes the 
distribution of wave sizes: t w = In (9/5)/ In 3. Numerically, it is found that the number of waves 
in an avalanche is proportional to the number of distinct sites toppled in the avalanche. This 
leads to a conjecture for the exponent r that determines the distribution of avalanche sizes: 
t = 1 + t w = In (27/5)/ In 3. Our predictions are in good agreement with the numerical results. 
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1 Introduction 



The sandpile model [ij plays a special role in non-equilibrium statistical mechanics since it is the 
standard example of a system that reaches a critical stationary state without the necessity of 
fine-tuning any parameters. This property, which was named self organised criticality, has been 
suggested to lie at the basis of the widespread occurence of fractals and l//-noise in natural systems 
P]. Though the sandpile model seems to have little relevance for the description of real sand, it 
has the advantage that it can be threated with analytical methods. Hence, it plays a role similar 
to that of the Ising model in the theory of equilibrium critical phenomena. 

In the ten years that have passed since its introduction the model has been intensively studied, 
mostly in two dimensions. Besides extensive numerical work ||, several theoretical developments 
have been made [4-8]. These have shown that many properties of the model are closely connected 
to those of spanning trees, which in turn can be described by the (q — ► 0)-limit of the Potts model 
H]. These connections have, for the case of two dimensions, led to exact results for the exponents 
that describe the distribution of sizes of waves and last waves || Q. Moreover, the relation to 
spanning trees has led Priezzhev and his coworkers || to a conjecture (again in d = 2) for the 
exponent governing for the sizes of avalanches, which in a way is the 'most wanted' of the sandpile 
exponents. The conjectured value is in good agreement with numerical estimates. In the present 
paper, we present results for the sandpile model on the Sierpinski gasket, a well known deterministic 
fractal. We have two main reasons for doing this. Firstly, we want to extend the theory for waves 
and the conjecture for avalanches to the case of the Sierpinski gasket and compare them with 
the results coming from extensive simulations. This can also lead to a better appreciation of the 
results in two dimensions. Secondly, the development of effective renormalisation group approaches 
is a major goal in the field of non-equilibrium statistical mechanics. In this respect deterministic 
fractals can play an important role since often real space renormalisation approaches can be worked 
out exactly for them. Some real space approaches to sandpile models have been introduced and the 
results presented here can be useful in testing their validity || [To| . We also mention that we are 
not the first to study the sandpile model on a Sierpinski gasket |LlJ . A previous study was however 
purely numerical and was moreover restricted to an investigation of avalanches. 

This paper is organised as followed. In section II we introduce the sandpile model, and give a brief 
overview of the existing theoretical work. In section III, we present our numerical results for waves 
and avalanches on the Sierpinski gasket. In section IV, our data are compared with exact results 
which we can obtain for wave exponents. Based on our numerical results, we also give a conjecture 
for the exponent r. Finally we present some concluding remarks in section V. 



2 Sandpiles, spanning trees and waves 



The (Abelian) sandpile model can be defined on any graph, but for definiteness we will introduce 
it in the contexts of the square lattice and of the Sierpinski gasket. In both cases each vertex of 
the lattice has four nearest neighbours. To each such vertex i we associate a height variable m 
which can take on any positive integer number. We also introduce a critical height z c of the height 
variable which we will take equal to four for all vertices. Furthermore, it is important to stress 



that we work on a finite part of the lattice (which we will denote by £), which contains N vertices. 
In the case of the Sierpinski gasket, N is trivially related to the number of iterations n used in 
constructing the fractal. 

The sandpile model has two time scales. On a slow time scale we drop sand at a randomly selected 
site and thereby increase the height variable by one unit: z% — > z% + 1. When at a given site, 
Zi > z c , we say that the site becomes unstable and perform a toppling of that site which involves 
the following operations: 



z-i 



Zi + A 



i-.i 



where 



A 



«-j 



^■j 



1 if i and j are nearest neighbours 
otherwise 



The matrix A is the discrete Laplacian. Through toppling, neighbouring sites can become unstable, 
topple themselves, create new unstable sites, and so on. This avalanche of topplings proceeds on 
the second, very fast, time scale, so that no new grains of sand are added before the avalanche is 
over. Sand can leave the system when a boundary site topples. Note that (see figure 1) in the 
Sierpinski gasket there are only three sites through which sand can flow out. An avalanche is over 
when all sites are stable again. We call also the resulting configuration of height variables stable. It 
is not difficult to see that the order in which unstable sites are toppled doesn't influence the stable 
configuration which is obtained when the avalanche is over. This Abelian nature of the sandpile 
model is crucial in the further analysis of its properties g . 

For further reference it also necessary to introduce the matrix G which is the inverse of A and 
which is therefore the lattice Green function of the Laplacian. It was shown by Dhar that the 
element Gij of the lattice Green function gives the expected number of topplings at site j after a 
grain of sand was dropped at site i. 



The criticality of the model manifests itself in power law distributions which are found for many 
quantities. For example, let D(s, N) denote the distribution of the number of distinct sites s toppled 
in an avalanche. It is found numerically that D(s,N) obeys a scaling law 

D(s,N) ^s~ T F{s/N) (2) 

Extensive simulations of the sandpile model in two dimensions lead to the estimate r = 1.22 Q. 
It is also found numerically that avalanches are compact, i.e. the fractal dimension D of the set 
of sites toppled in a avalanche is equal to that of the lattice on which the model is studied. In 
particular, we have for the case of the Sierpinski gasket of figure 1, D = In 3/ In 2. 

The number of stable configurations of the model equals z^ . However, it has been shown by Dhar 
Q that the number of recurrent configurations in the stationary state is exponentially smaller. In 
fact, there is a one-to-one correspondance between these recurrent configurations and the number 
<Sjv of spanning trees on a lattice £* which consists of C extended with an extra site, known as the 
sink, which represents the boundary j5|. It is known that Sn grows in general as ji N+ °^ N \ where \i 



is called the connective constant of spanning trees |[^]. Spanning trees are described by the g-state 
Potts model in the limit q — » [13]. The same model also describes resistor networks. 



Because of the Abelian nature of the sandpile model, we can perform the toppling of unstable sites 
in any desired order. A particularly interesting order leads to the concept of waves j^. Suppose 
that an avalanche starts at a site iq and that after some topplings iq becomes unstable again. One 
can then keep io fixed for the moment and continue with the toppling of other sites untill all of 
them (apart from iq) are stable. It is easy to show that in such a sequence all sites topple at most 
once. This set of topplings is called the first wave. Next, we topple the site i$ for the second 
time. If after some topplings it becomes unstable again, we keep it fixed, and topple all the other 
unstable sites. This set of topplings constitutes the second wave. We continue in this way untill we 
finally reach a stable configuration. We can in this way decompose any avalanche in a set of waves. 
Of course, subsequent waves in an avalanche are highly correlated, but we can also consider the 
statistical properties of the set of all waves, irrespective of the avalanche in which they occur. In 
this way, we can look at the distribution D(s w ,N) of the number of sites s w toppled in any wave. 
Again, it is found that D(s w ,N) has a scaling form 

D(s w ,N)^s- T -F(s w /N) (3) 

where we have introduced the exponent t w for waves. 

The properties of the last wave in a given avalanche are also of interest. The distribution of the 
number of sites si w toppled in this last wave leads to the introduction of still another exponent, 
denoted as r\ w . Because s\ w < s, we arrive at the obvious upper bound 

r < t 1w (4) 

It turns out that the properties of waves are easier to analyse than those of the avalanches. Without 
going into details we summarize the following important results: 

• There is a one-to-one correspondance between waves and two-rooted spanning trees (one root 
corresponds with the sink, the other with the site io) on C* |J7|]. Furthermore, the element Gij 
of the lattice Green function is given by the ratio of the number of two rooted spanning trees 
(in which i and j are in the same subtree) to the number of (one rooted) spanning trees Sjy. 
Hence, properties of waves can be determined from G. From the known asymptotic behaviour 
of the lattice Green function on an Euclidean lattice, Gij ~ \i — j\ 2 ~ , one then finds 

2D -2 

Tw = D (5) 

• The exponent t\ w can be obtained from the probability that b sites get disconnected from a 
spanning tree when a randomly chosen bond of the tree is deleted. In particular, one finds 
that T\ w can be expressed in terms of the chemical dimension z of spanning trees || . On an 
Euclidean lattice the precise result is 

r lw = l+ 2 -^- (6) 



From these results, one obtains, in D = 2 (where z = 5/4 Q) 

t w = 1 r lw = 11/8 

A final theoretical development which we need to discuss leads to a conjecture on the exponent r 
||]. Priezzhev et al. first relate r to another exponent, a. To introduce this exponent, it is assumed 
that the number of sites toppled in consecutive waves (say the fc-th and k + 1-th wave) decreases. 
The change in size, As w = s w (k) — s w (k + 1) is assumed to occur in a self-similar way, so that 

As w ~ si (7) 

One can then immediately show || that 

r + a = 2 (8) 

Finally, using the relation between waves and two-rooted spanning trees, and considering possible 
processes by which the size of a wave can decrease, it is found that a = 3/4 and hence r = 5/4. 
This prediction is in excellent agreement with numerical estimates, and may very well be exact. 



3 Numerical results 



In order to obtain numerical estimates for the exponents r, t w , T\ w and a on the Sierpinski gasket 
we have performed extensive simulations of the sandpile model on this fractal structure. In these 
simulations, we have considered lattices with n = 3 up to n = 8 (the latter contains N = 9843 
sites). For each of these we have studied more then (1000 N) avalanches. 

In figure 2, we show our results for the distribution of the size of the waves. As can be seen, log- 
periodic oscillations are superposed on the power law decay of the distribution. The behaviour of 
D(s, N) can therefore be described by complex critical exponents. This is a well known phenomenon 
in systems with discrete scale invariance JL4|]. The imaginary part of the critical exponent is in a 
sense trivial since it can be related to the discrete rescaling factor of the fractal lattice. From a 
careful analysis of the data one can hence obtain finite size estimates of the exponent T w {n) as a 
function of n. These can be extrapolated by making the usual assumption 

C 

Tw (n) =t w + (9) 

inn 

In figure 3 we show this extrapolation. It leads to the estimate 

t w = 0.47 ± 0.05 (10) 

It is difficult to make a precise determination of the error in this estimate. The error we give is 
therefore more indicative. 

Figure 4 shows our results for the distribution of the size of the last wave. From an examination 
of these data, we find 

t 1w = 1.66 ± 0.05 (11) 



Finally, our data for the distribution of avalanche sizes, shown in figure 5, yield the exponent 
estimate 

t = 1.46 ±0.05 (12) 

in agreement with the result found by Kutnjak-Urbanc et al. j[l|], who obtained r = 1.51 ± 0.04. 

We have also investigated the behaviour of As w in order to determine the exponent a. To obtain 
this exponent, we have to consider only those cases in which As w > (in fact, there are many cases 
for which the opposite inequality holds). In figure 6 we show our data (after integration over bins). 
As can be seen from the figure, it is not easy to determine a clear value of a. Our best estimate is 

a = -0.15 ±0.1 (13) 



In the next section, we will give a lower bound for a. Our estimate (13) is very close to that bound. 
It is however difficult to imagine how a can be negative. We therefore turned to an alternative way 
to determine a. This is done by investigating how the average number n w of waves in an avalanche 
scales with the size of the avalanche. When the size of consecutive waves decreases, the size of the 
first wave is proportional to the size of the whole avalanche. From (|7|) we then obtain 

n w ~ s 1 "* (14) 

In figure 7 we show data for this quantity for n = 6 and n = 8. 

We find very clear evidence for a linear relation between n w and s. This relation breaks down only 
for the largest avalanches, which are those that span the lattice. But, as pointed out by Kutnjac- 
Urbanc et al. [11], these avalanches are somewhat peculiar. They occur here because only through 
such avalanches is it possible for sandgrains to leave the system. 

From the linear relation we find 

a = 0.0 ±0.01 (15) 

These data lead us to conjecture that a = exactly. As we will see in the next section, this result 
will be very useful in determining the exponent r. 



4 Exact and conjectured results 



We now present some exact results for the connective constant /i for spanning trees ( and hence for 
the number of recurrent states) and for the exponents t w and t\ w . We make also a conjecture for 

T. 

As discussed in section 2, there exists a close connection between the theory of the sandpile model 
and the (q — > 0)-limit of the Potts model. This connection should allow the determination of several 
properties of the sandpile from a study of that model on the Sierpinski gasket. For example, the 
connective constant of the spanning trees can be obtained from the free energy of the Potts model 
when q — > 0. This free energy in turn can most easily be obtained by performing a simple but 
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exact real space RG calculation for the Potts model on the Sierpinski gasket. Consider therefore 
the g-state Potts model for which we will denote the spin variable at vertex i by Ui. On a Sierpinski 
gasket, these spins can also be labeled according to the 'generation' in which they appear during 
the iterative construction of the fractal. The reduced Hamiltonian H(K) of the Potts model is 
given by 

H(K) = KJ2^ (16) 

(id) 

A renormalisation group for H(K) on the Sierpinski gasket can be obtained by performing a partial 
trace, denoted by Tr* over the spin variables which have the same 'generation label': 

e H(K') + g(K) =Tr * e H(K) (17) 

The free energy of the model can then be obtained by a summation of g(K) along the renormali- 
sation group trajectory. We have performed such a calculation. After sending q — ► the RG flow 
is remarkably simple 

K' = -K (18) 

5 



MK) = 50 ^3 g 3/2 (19) 



from which we obtain the connective constant for spanning trees after a straightforward calculation. 
The result is 

,1/3 



/i = 50 1/d A/3/5« 2.854 (20) 



As mentionned in section 2, the exponent n w can be obtained from the probability P s {b) that upon 
deleting a bond at random from a spanning tree, a subtree of size b gets disconnected. On an 
Euclidean lattice this relation leads to the result @. On a fractal lattice, the relation has to be 
modified. The correct expression now also involves the fractal dimension D w of a random walk on 
the fractal. This was shown in reference [ffBIl, where it was found that 



P s (b) ~b- Dw/z 
Then, repeating the arguments in we obtain 

ri w = 1 + ^p (21) 

This is a quite general result for ti w which should for all cases (fractal and Euclidean). 

For spanning trees on the Sierpinski gasket, z was recently determined by Dhar and Dhar in a 
study of loop-erased random walks [fi^] . Their result is 

z = ln[(20 + V205)/15)] 

In 2 { } 

[Before proceeding, we would like to point out a general recipe for determining z, using again an 
RG approach. For a tree, the chemical dimension equals the red bond dimension Dr. The red 



bonds of a connected graph are those edges which when cut, disconnect the graph. Coniglio [16 1 , 
and independently Saleur and Duplantier |l7j], pointed out that the red bond dimension for the 
critical Potts model can be obtained by adding a term 

hi 

to the Potts hamiltonian (|l6|). The red bond dimension can then be obtained from the linearised 
RG equation for M at the critical Potts fixed point. Hence, it is possible to get z for spanning trees 
by performing an RG calculation (either exact or approximate) on an extended Potts model. We 
will present results from such an approach in a future publication.] 

Inserting ( f22|) in ( |2l| ) and using the value of D w for the Sierpinski gasket (D w = In 5/ In 2, |]l^]) we 
get 

In [(20 + ^205)775] 
Ino 

This prediction is in good agreement with the numerical value found in section 3. 

Next, we turn to the exponent t w . This should be determined from the behaviour of the Green 
function. We will assume that the Green function G{j decays for large r = \i — j\ as 

Gij ~ r~ 2x (24) 

The Green function gives the expected number of topplings at site j when a grain of sand has been 
dropped at site i. Hence it consists of all waves that cover the distance between site i and j. Since 
waves are compact, we obtain 



/•oo 

Gij ~ / D(y)dy 



where we have neglected finite size effects. From fl24|) we then get a general expression for the 
exponent t w 

2x 

t w = 1 + — (25) 

To determine t w for the particular case of the Sierpinski gasket, we only need to determine the 
value of x. This exponent was determined in early studies on the behaviour of the resistance on 
fractals with the result [|l9| ] 

2x ^ln(5/3) 



This leads to the prediction 



In 2 



n (9/5 . . 

; ' ' w 0.535 (26) 

In 3 



This value is again in agreement with the one determined numerically. In principle, the exponent 
x can also be obtained from the q — > limit of the Potts model. In fact, it can be obtained from 
the renormalisation group equation for K, see (|X 
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Finally, we turn to the exponents r and a. It is easy to generalise the result @ which was obtained 
in D = 2. The size s w of any wave must be less then or equal to the size s of the avalanche in which 
it occurs. Moreover, the number of waves whose size lies between s w and s w + ds w is according to 
(|7|) proportional to s~ a . Combining these two results we obtain for the size distribution of waves 

DM ~ s- a - T - 1 

from which we immediately obtain the appropriate generalisation of (||) 

a + t = 1 + t w (27) 



Our numerical results (|To|), (|T^) and (15) are nicely consistent with this picture. From the upper 
bound (H) , we obtain a bound on a 

a > 1 + t w - T iw ss -.177 

Our first estimate of a, given in (|13|) , is just consistent with this bound. 

Most interestingly our conjecture that a = leads to a conjecture for r 

In (27/5) , x 

t = \ + t w = V /; ^1.535 28 

In (3) 

which is in excellent agreement with our numerical results. 

We have also tried to extend the argument of Priezzhev and his coworkers for the exponent a to 
the fractal case. As far as we could see, such an argument always leads to a value of the exponent a 
which is of order one, and which therefore is completely in disagreement with the numerical result. 
It therefore seems that the way in which consecutive waves are related is different for the fractal 
case than for the Euclidean one. 



5 Discussion and conclusions 



In this paper, we have presented extensive simulations for waves and avalanches on the Sierpinski 
gasket. One of the most surprising results is the fact that the exponent a is very close to zero. 
This implies that the way in which waves become smaller is almost independent of their size. This 
is most probably a direct consequence of the fractal nature of the Sierpinski lattice. 

We have also extended the existing theory for waves to the present case. This has led to exact 
expressions for the exponents ri w and t w which are close to the numerically determined values. 



Most interestingly, assuming that a is zero, we are led to the conjecture (28) for the exponent r. 
Also this prediction is in good agreement with our simulation results. 

In this way, there is now a second situation (after the two-dimensional case) for which wave expo- 
nents are determined exactly and for which a conjecture for the r-exponent exists. 
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The result a = indicates, in our opinion, that there is good hope that the sandpile model on a 
Sierpinski gasket may be solved exactly. We are currently using real space RG approaches which 
have been introduced for sandpile models ||, [jl| , to see whether they may be worked out exactly 
on this lattice. 
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Figure 1: The simple Sierpinski gasket with n 
three sites indicated by an arrow. 



2. Sand can only leave the lattice through the 
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Figure 2: Distribution of wave sizes. The data give numerical results for n = 3 to n = 8 (top to 
bottom). 
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Figure 3: Extrapolation of T w (n) (L = 2 n ). 
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Figure 4: Distribution of last wave sizes. Data are given for n = 3 to n 



15 



0.1 
D(s) 

0.01 



0.001 



0.0001 



10 



-5 



n i mm 



<2g?. ~. V^T*j 



^* 



VVtfV \j* £ 



J I I I III 




10 



100 



1000 s 10 4 



Figure 5: Distribution of avalanche sizes. Data are given for n = 3 to n 
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Figure 6: Decrease in size (averaged over bins) of consecutive waves, versus wave size, for n = 6 
and n = 8 (upper dataset). 
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Figure 7: Number of waves in an avalanche versus size of the wave, for n = 6 and n = 8 (upper 
dataset). 
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